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ABSTRACT 

We analyze A^-body simulations of halo mergers to investigate the mechanisms responsible for driving mix- 
ing in phase-space and the evolution to dynamical equilibrium. We focus on mixing in energy and angular 
momentum and show that mixing occurs in step-like fashion following pericenter passages of the halos. This 
makes mixing during a merger unlike other well known mixing processes such as phase mixing and chaotic 
mixing whose rates scale with local dynamical time. We conclude that the mixing process that drives the sys- 
tem to equilibrium is primarily a response to energy and angular momentum redistribution that occurs due to 
impulsive tidal shocking and dynamical friction rather than a result of chaotic mixing in a changing potential. 
We also analyze the merger remnants to determine the degree of mixing at various radii by monitoring changes 
in radius, energy and angular momentum of particles. We confirm previous findings that show that the majority 
of particles retain strong memory of their original kinetic energies and angular momenta but do experience 
changes in their potential energies owing to the tidal shocks they experience during pericenter passages. Fi- 
nally, we show that a significant fraction of mass 40%) in the merger remnant lies outside its formal virial 
radius and that this matter is ejected roughly uniformly from all radii outside the inner regions. This highlights 
the fact that mass, in its standard virial definition, is not additive in mergers. We discuss the implications of 
these results for our understanding of relaxation in collisionless dynamical systems. 

Subject headings: cosmology: theory — dark matter:halos — galaxies: relaxation — halos: structure — 
methods: numerical 



1. INTRODUCTION 

Almost four decades have passed since "violent relaxation" 
was first proposed by Lvnden-Bell ( 1967) as a mechanism to 
explain the remarkable regularity and smoothness of the light 
distributions in elliptical galaxies. Several au thors had noted 
' that t he time scale for two-body relaxation dChandrasekhaii 
■ [1941 was several orders of magnitude too long to be rele- 
] vant for the evolution of galaxies (Zwicky 1939). Lynden- 
Bell derived the distribution function that would result from 
rapid fluctuations in the potential of the galaxy, such as those 
] that might occur during gravitational collapse or following a 
merger. The process of violent relaxation outlined by Lynden- 
Bell produces a smooth distribution function as a result of 
the gravitational scattering of particles by the time depen- 
dent global gravitational potential. He argued that this process 
could result in a mass-independent and smooth final distribu- 
tion of particles in phase-space. 

Ano ther process that was considered by iLvnden-Belll 
(11967) is "phase mixing" - a process by which an initially 
compact ensemble of non-interacting phase-space points in a 
time-independent potential can be stretched into a long ribbon 
in phase-space as a result of small differences in the initial 
energies (or angular momenta) of particles in the ensemble. 
Phase mixing conserves the fine grained distribution function 
but in a coarse grained sense it produces smooth/relaxed look- 
ing distributions on timescales that depend only on the range 
of values of the orbital integrals (or actions) in the ensemble 
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dBinnev & Tremainelll987l) . 

In the last decade, a new type of mixing process, which 
occurs due to the presence of large fractions of chaotic or- 
bits in a galactic potential (such as when box orbits are 
scattered by a central supermassi ve black hole in a triax- 
ial galaxy), has b een inv estigated (K andrup & Mahon 199^ 
Merritt & Vallu ril 119961: iHabib et al. .1997; Kandrup et al] 
2000t [Kandrup & SiopisI 120031) . These studies have shown 
that microscopic ensembles of iso-energetic test particles 
which occupy an infinitesimal volume of phase-space sur- 
rounding a chaotic orbit can evolve rapidly to a near-invariant 
distribution that uniformly covers the energy surface available 
to the ensemble. This process of diff'usion to a near-invariant 
distribution is termed "chaotic mixing" and has been shown 
to occur in a variety of galaxy-type potentials. These authors 
argued that in a triaxial galaxy with a supermassive central 
black hole, a significant fraction of orbits are chaotic, and con- 
sequently, this mixing process could drive secular evolution to 
more axially symmetric distributions. 

The physically interesting aspect of this mixing process is 
that the evolution to a near-invariant distribution has been 
shown to occur on timescales of order 5 to 100 orbital crossing 
times in strongly chaotic systems making it a potentially im- 
portant pr ocess for the evol ution of galaxies. This was demon- 
strated by iMerritt & Ouin lan (1998) who employed A^-body 
simulations of triaxial elliptical galaxies containing growing 
supermassive central black holes and found the triaxial figures 
to evolve to more axisymmetric shapes. 

Like phase mixing, chaotic mixing conserves the fine- 
grained phase-space density. Unlike phase mixing, its rate 
does not depend on the initial spread in the energy distribu- 
tion of particles in the ensemb le, but on how strongly chaotic 
the orbits in the ensemble are (IMerritt &Vallurilll996l) . 

However two concerns remain regarding these earlier ex- 
periments. First the rate of evolution to the invariant measure 
depends on the properties of the ensemble of particles. For 
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ensembles near a "sticky" or weakly chaotic orbit, a near- 
invariant di stribution may not be a ttained in the lifetime of 
the galaxy jMerritt & Valluril 11996 ). Second, the majority 
of the numerical studies which found strong chaotic mix- 
ing and short time scales for evolution were not carried out 
in self-consistent potentials leading to concerns that in self- 
consistent systems such mixing could be self-limiting. 

Since all orbits in a time-dependent potential with long- 
lived oscillations can, in principle, exchange energy with 
the potential, a large fraction of the orbits will not con- 
serve the energy integral. In such systems previous authors 
have found that a sig nificant fraction of orbits can be chaotic 
(iKandrup et alJl2003l) . This raises the possibility that the un- 
derlying physical process driving relaxation in violently fluc- 
tuating potentials could be a consequence of chaotic mix- 
in g. The first step in m aki ng this connect ion was studied 
bv lKandrup et all J2003l) and iTerzic & Kandr up (2004), who 
demonstrated that large fractions of chaotic orbits arise in po- 
tentials subjected to long-duration (damped) periodic oscilla- 
tions. These authors found that the fraction of chaotic orbits 
depended on the frequency as well as the amplitude of the os- 
cillations in the potential and attributed the chaos to a "broad 
parametric resonance". 

In this paper we analyze the mechanisms driving the ap- 
proach to equilibrium following the collisionless merger of 
two self-gravitating systems. Our principal interest is in un- 
derstanding whether changes in the fine-grained distribution 
function (e.g., changes in the angular momenta and energies 
of nearby particles) occur as a result of exponential instabili- 
ties of orbits or some more straightforward energy exchange 
process which transfers kinetic energy of the merging halos to 
the internal energies of particles. 

We consider two dynamical heating effects that are not usu- 
ally discussed in this context. The first is the effect of com- 
pressive tides that arise when one collisionless gravitating 
system passes through another on a time scale that is short 
compared to the internal dynamical time of the infalling sys- 
tem. Such tides can impulsively heat particles as a result of 
the transient deepening of the net potential. Impulsive com- 
pressive tidal shocks have been studied previously and are 
well known to produce large chan ges in the internal structure 
and evolution of globular clusters (ISpitzer & ChevalieJll973l : 
iGnedin et al."1999.). galaxies orbiting within a cluster poten- 
tial d^lluri 199 3b, and DM subhalos orbiting with in larger 
hosts lKravtsov et al.ll2004l: iKazantzidis et al J'2004bV 

The second process is dynamical friction ( Chandra sekhaii 
Il943h . which is the primary mechanism that causes the forma- 
tion of a single remnant when two collisionless halos merge. 
When collisionless halos merge, dynamical friction transfers 
their relative linear and angular orbital momenta to the indi- 
vidual particles. The original version of the Chandrasekhar 
( Il943h dynamical friction theory applied to the deceleration 
of a massive star moving in an infinite homogenous isotropic 
stellar distribution. However several studies have shown that 
with minor modifications the theory is applicable in systems 
that resemble real galaxies and clusters (White 1976; Binnev 
ll977H Velazauez & White 1999). Dynamical friction is found 
to be enhanced in systems with inhomogenou s density pro- 
files (such as cuspy profiles iDel Popololl2003h . as compared 
with systems with uniform homogenous distributions. Recent 
A^-body studies of live satellites in exte nded A^-body galaxies 
(|Jiang & Binney 2000; Fujii et al. 2006) have shown that the 
timescale for dynamical friction is shorter than that predicted 
by Chandrasekhar due to two additional components of drag 



force: the tidal torque from the leading tidal tail and an en- 
hanced wake due to particles that are no longer bound to the 
satellite but still trail behind it in the same orbit. 

It is important to emphasize at the outset, that the terms 
"mixing" and "relaxation" are often used rather loosely (and 
often interchangeably) in the astrophysics literature to mean 
the approach of a gravitating system to global as well as small- 
scale equilibrium. The most restrictive definition of "relax- 
ation" is a process which leads to the complete "loss of mem- 
ory of initial conditions" or loss of correlations between ini- 
tially nearby particles in phase space. The large changes in 
integrals of motion required for such relaxation only occur 
as a result of two-body interactions (IChandrasekhad [T943h 
and as a result of "violent relaxation". Both processes cause 
changes in the fine-grained distribution function. In contrast 
phase mixing and chaotic mixing in time-independent poten- 
tials, under the right conditions, can lead to large changes in 
correlations but conserve integrals of motions (if they exist) 
and conserve the fine-grained distribution function. In this 
paper we focus our attention on loss of correlations in phase 
space accompanying a collisionless merger and use the term 
"mixing" to connote the approach to macroscopic equilibrium 
and the accompanying change in the fine-grained distribution 
function (even though this is quite small). 

According to the prevailing cosmological model for struc- 
ture formation, galaxies and dark matter (DM) halos in the 
Universe form hierarchically via multiple mergers. Numerous 
minor mergers and several major mergers of subhalos eventu- 
ally form the dark halos at 2; = 0. Thus, the merger history 
of cosmological structures is complex, and DM halos are pre- 
dicted to abound in substructure and have mass profiles that 
may still be evolving today. A detailed analysis of the physics 
of gravitational relaxation, with the specific goal of under- 
standing what processes drive the DM and stellar component 
to equilibrium, can be more readily done using dissipationless 
controlled merger simulations. 

Several aspects of the universality in the structure and 
properties of dark matter halos have lead to new interest in 
the detailed physical processes driving gravitational relax- 
ation. The discovery of "universal" dark matter halo pro- 
files resulting from c osmolo gical A^-body simulations (e.g., 
iDubinski & Carlberd [19911: iNavarro et ail [19961 \l99% . the 
more recent realization that these profiles have only a "nearly 
universal" form with a correlation betwe en the density pro- 
file's shape parameter and the halo mass dMerritt et al.ll2005l 
I2OO6I) as well as the finding that phase-space density profiles 
of these halos have power-law distri butions over about 2.5 or- 
ders of magnitude in radial extent (iTavlor & Navarroll2001l : 
lArad et al.l2004HAscasibar & Binnevl2005h all point to a uni- 
versality in the processes driving structure formation. In par- 
ticular, violent relaxation and the statistical distribution func- 
tions th at result from this process are being actively inves- 
tigated dXrad & Lvnden-Bell 200l lArad & Johanssonll2005l: 
iDehnenI 2005h . Several studies of the remnant density pro- 
files have shown that the remnants retain a surprisingly strong 
memory of their initial density profiles despite the strength of 
the relaxati on proces s (Barnes [ 19961: iBoylan-Kolchin & Mai 
.2004 ; Kazantzidis et al. 2 006 ) . 

Several authors have criticized Lynden-Bell's 1 19671 origi- 
nal suggestion that it is potential fluctuations that occur dur- 
ing gravitational collapse that drive the system to equilib- 
rium, on the grounds that time dependence of the potential 
alone merely causes particles to be relabeled in energy and 
causes no mixing of nearby particles in energy and angu- 
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lar momentum. It has been therefore argued that a process 
like chaotic mixing is actually responsible (at the microscopic 
level) for d riving the evolutio n of the s ystem to equilibrium 
dMen-itt & Valluri 1996; KandruE.etafl^003). 

An analysis of the merging process can help us better un- 
derstand the physical mechanisms that drive the evolution to 
equilibrium and cause the redistribution of nearby particles 
in phase-space. In this investigation we restrict ourselves to 
the study of mixing in positions, velocities, energies and an- 
gular momenta of particles. The goals of this study are: (a) 
to determine whether the relaxation of a remnant to an equi- 
librium distribution following a coUisionless A^-body merger 
is driven by chaotic mixing that is believed to take place in 
time-independent potentials (b) to quantify the degree of re- 
laxation/mixing in various phase-space parameters that occurs 
during the merger by examining the remnants of major merg- 
ers of DM halos, and (c) to obtain a better understanding of 
why steep cusps are as robust as they are. Both the time- 
dependent nature of the problem, as well as the fact that A^- 
body systems exhibit exponential sensitivity to initial condi- 
tions (as discussed in greater detail in § I2.2| |. make the first of 
the above goals particularly challenging. 

The paper is organized as follows. In § 2 we describe the 
A^-body experiments that were conducted and define several 
new measures that we use to quantify the rate of separation 
of nearby orbits and the degree of "mixing" or "relaxation" 
in phase-space. In § 3 we discuss the results of the applica- 
tion of our measures of mixing to the merger experiments by 
considering the separation in various phase-space variables of 
pairs of nearest neighbors in phase-space (§ 13. lb as well as by 
monitoring the rate of mixing of ensembles of nearby parti- 
cles (§ I3.2| |. In § 4 we examine how particles are redistributed 
in phase-space following the merger. In § 5 we summarize 
our results and discuss the implications of our findings for 
understanding relaxation during mergers and the mixing that 
it produces. 

2. NUMERICAL METHODS 

We investigate the mechanisms responsible for violent re- 
laxation by analyzing A^-body simulations of binary mergers 
between dark matter halos. The halo models follow density 
profiles that are describe d by the g eneral (a,/?, y) spherical 
density law (eq 1.) (e.g., Hernquist 1990; Zhao 1996), where 
y denotes the asymptotic inner slope of the profile, fi corre- 
sponds to the outer slope, and a determines the sharpness of 
the transition between the inner and outer profile. 

We consider two main halo models specified by particular 
choices of the parameters a, y6, and j. The first model fol- 
lows the Navarro et al. (1996, hereafter NEW) profile (with 
( a, yS, y) = ( 1 , 3 , 1 )), while the second model (with ( oc, j0, y) - 
(2,3,0.2)) corresponds to a profile with a shallower inner 
slope. A^-body halo models are constructed using the exact 
phase-space distribution function under the assumptions of 
sph erical symmetry and an isotropic velocity dispersion ten- 
sor dKazantzidis et al.ll2004al) . Each of the initial DM halos 
has a virial mass of Mvir = lO'^M© implying a virial radius of 
Tvir - 256.7 kpc, and a concentration of c = 12, resulting in a 
scale radius of ^ 21 .4 kpc. It is worth emphasizing that the 
adopted value of Mvir serves merely practical purposes and 
does not imply anything special about the particular choice of 
mass scale. Hence, our conclusions can be readily extended 
to mergers between equal-mass systems of any mass scale. 



We analyze three different merger experiments which were 
conducted w ith the multi- stepping, parallel, tree A^-body code 
PKDGRAV (IStade]||2001h . PKDGRAV uses a spHne soften- 
ing length, such that the force is completely Keplerian at twice 
the quoted softening length, and multi-stepping based on the 
local acceleration of particles. The first experiment followed 
the encounter of two NEW halos (referred to hereafter as run 
Bpl), while in the second experiment an NEW halo merged 
with a halo having an inner slope of y = 0.2 (referred to here- 
after as run hBpl). The initial halo models in these two runs 
were sampled with A^ - 2 x 10^ particles and forces were 
softened with a spline gravitational softening length equal to 
e = 1 .5 kpc. In order to minimize any concern that our results 
might be compromised by numerical resolution effects, we 
analyzed an additional merger experiment between two NEW 
halos increasing the mass resolution by a factor of 10 and scal- 
ing down the softening lengths according to e oc A^^'^^ (here- 
after referred to as run HRBp). Initial conditions for binary 
mergers were generated by building pairs of halo models and 
placing them at a distance equal to twice their virial radii. In 
this study we only discuss mergers of systems on parabolic or- 
bits owing to the fact that this particular orbital configuration 
is the most t ypical of merging h alos in cosmological simu- 
lations (e.g.. IZentner et all 2005h . These merger simulations 
(labeled Bpl, hBpl, and HRBp) were analyzed and described 
in greater detail in Kazantz idis et al. (2006 ). Exten sive con- 
vergence tests carried out bv lKazantzidis et al.l ( l2006h indicate 
that isolated halo models did not deviate from their original 
distribution functions on timescales as long as 100 crossing 
times even at small radii (r ~ e). 

In order to distinguish the mixing that results from the ex- 
ponential instability of the A^-body problem (see § I2.2l i from 
the effects of the mixing resulting during the merger, we also 
performed A^-body simulations of the evolution of a spherical 
isolated NEW halo. All orbits in an isolated spherical halo are 
expected (in the smooth potential or large A^ limit) to be "reg- 
ular orbits" that conserve four isolating integrals of motion. 
This simulation therefore serves as an important control ex- 
periment with which to compare the evolution in the merger 
simulations. 

2.1. Macroscopic Evolution of Merger Remnants 

A self-gravitating distribution of particles that is out of 
equilibrium will experience potential fluctuations, and the po- 
tential and kinetic energies of particles will oscillate back and 
forth. This behavior is governed by the time-dependent virial 
theorem: 

,2, 



-/ 

2 



2r + y 



where / is the moment of inertia, T is the total kinetic en- 
ergy, and y is the total potential energy of the system of par- 
ticles, with all quantities defined with respect to the center of 
mass of the system. Eor a time-independent gravitating sys- 
tem, / = = 2r + y. A robust quantitative measure of how far 
a gravitating system is from equilibrium can be obtained from 
the virial ratio 2r/|y|, which is unity for a system in global 
dynamical equilibrium. 

In Eigure[T] the top panel shows the separation of the most- 
bound-particle (MBP) of each halo in the merger, the middle 
panel shows the evolution of the virial ratio 2T/\V\ and the 
lower panel shows the change in total potential energy of the 
system {V), as a function of time. We note that pericenter 
passages (seen as minima in the MBP separation) correspond 
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Fig. 1. — Top: The separation between the most-bound particles of the two 
merging NFW halos (run Bpl) as a function of time; Middle: time evolution 
of the virial ratio 2T/\V\: Bottom: time evolution of the total potential energy 
of the system V. 

to maxima in the virial ratio and minimum in potential energy. 
For f < 7 Gyr the virial ratio 2T/\V\ and potential V undergo 
strong fluctuations but remains close to unity thereafter We 
refer to this time t ~ 7 Gyr to be the time when the system is 
globally relaxed. 

2.2. Microscopic Chaos or the Miller Instability 

In real galaxies, there are two components of the gravita- 
tional accelerations on particles, a component arising from the 
background potential which is smoothly varying with time, 
and a second component arising from the discrete nature of 
stars (and possibly dark-matter particles). The importance of 
the discret e component is defined by the two-body relaxation 
timescale (iChandrasekhad 1 1 943t iBinney & Tremaind fT987h 
and is larger than the age of the Universe in galaxy-sized sys- 
tems. Therefore, it is generally assumed that equilibrium and 
non-equilibrium evolution of the phase-space density distribu- 
tions in galaxy-size objects can be studied in the mean-field 
limit, where the collisionless Boltzmann equation can be ap- 
plied. 

The Lyapunov exponent (iLichten berg & Liebermanll 19921) 
is one of the most popular ways of measuring the chaoticity 
of an orbit in a smooth potential. It measures the rate of diver- 
gence of two infinitesimally nearby orbits in phase-space. If 
the rate of divergence of phase-space points is linear, then the 
Lyapunov exponent A = and the orbit is termed regular. In 
a smooth equilibrium potential with three-spatial dimensions 
such orbits will conserve at least three isolating integrals of 
motion. If, on the other hand, the rate of divergence of the 
phase-space points is exponential, the orbit is termed chaotic. 



Thus, the Lyapunov exponent defines the time scale on which 
infinitesimally nearby trajectories in phase-space undergo ex- 
ponential separations in their phase-space coordinates. 

However, as was first noted by Milled d 19641) . the A^-body 
problem is chaotic in the sense that the trajectory of the 6A^- 
dimensional phase-space coordinate of the system exhibits 
exponential sensitivity toward small changes in initial condi- 
tions. This exponential instability (referred to in the literature 
as the "MiUer in stability"; e.g., He msendorf ^Merrittil2002l 
or "microchaos" ISideris & Kandrupll200 2) has been investi- 
gate d extensively i n several studies over the past three decades 
(see lMerrittll2005l for a review). These studies have shown 
that the largest A^-body Lyapunov exp onent A does not con- 
verge toward zero as increases ( Kandrup & Smiflj[T99Tal: 
iGoodman etanil993h. even for A^-body re alizations of inte- 
grable potentials (iKandrup & Siderisll2001h. 

In fact more recently iHemsendorf & MerrittI (j2002) have 
shown that the direct A^-body problem (with A^ < 10^) is in- 
herently chaotic and that the degree of chaos, as measured by 
the rate of divergence of nearby trajectories (with the Lya- 
punov exponent A), increases with increasing A^ with charac- 
teristic e-folding time equal to ~ 1/20 of the system cross- 
ing time (in the absence of softening). It has also been 
shown that the mean (e-folding time) increases with in- 
creasing particle number in the case of softened potentials 
jEl-Zant 2002) so that discreteness effects are reduced as 
the softening becomes com parable to interparticle separations 
(IKandrup & Smith|[T99Tbl) . 

The exponential separation of particles in energy occurs on 
the Miller instability timescale (f^,) which is much smaller 
than the two body r elaxation time (f,) (Kandrup et al. 1994| 
iHut & Heggijl200ll) . The diver gence of particles in energy 
saturates after a few t^ and then varies with time on the stan- 
dard two-body relaxation timescale indicating that the Miller 
instability is not an important process in changing the fine- 
grained distribution functions of gravitating systems. 

Despite the fact that particles separate exponentially in 
phase space, as A^ increases orbits in these systems become 
more and more regular, acquiring orbital characteristics closer 
and closer to those in smooth potentials, and the macro- 
scopic sep arations of orbit s satura te at smaller and smaller 
distanc es dValluri & Merritti I200a IKandrup & SiderisI 120011: 
ISiderisllIOOir Additionally, decades of A^-body integrations 
have demonstrated that, in many ways, the behavior of large- 
A^ systems matches expectations derive d from the co Uision- 
less Boltzmann equation (e.g., Aarseth & Lecajri975h . 

The fact that in an A^-body system the Lyapunov exponent 
will primarily measure the "Miller instability" and not macro- 
scopic chaos (arising from the background potential) means 
that it is has limited usefulness for measuring chaotic be- 
havior in A^-body systems. Other measures of chaos, such 
as those that measure changes in orbital characteristics (for 
instance, fundamental oscillation frequencies of orbits ; e.g., 
lLaskari[l99(illLaskar et alJlT99l I Valluri & Merrittl[T998h . rely 
on the ability to identify such frequencies from oscillations 
that last between 30-100 orbital periods. Recently, a new 
technique based on pattern recognition of orbital character- 
istics has been developed and is applicable in time-de pendent 
system s whose oscillations last 10-30 crossing times dSiderisI 
|2006|) . To our knowledge, none of the standard measures of 
chaos are applicable to A^-body orbits in potentials with non- 
periodic potential fluctuations lasting only a few (< 10) cross- 
ing times. This makes it necessary to define new ways of 
quantifying chaos and mixing. 
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2.3. Definitions of mixing in N-body systems 

In this paper we estimate the degree of mixing at various 
stages of the evolution of the merging systems (a) by consid- 
ering pairs of nearby particles in the A^-body simulation and 
tracing their separation in four different phase-space quanti- 
ties and (b) by carrying out mixing experiments that moni- 
tor the rate at which ensembles of 1000s of nearby particles 
evolve with time and reach a "near invariant" distribution in 
energy and angular momentum. 

2.3.1. Separation of pairs of nearby particles 

Several previous studies of this type have focused on the 
configuration space separation Ar of a pair of infinitesimally 
nearby test particles in a frozen A^-body potential. These 
studies picked pairs of non-self-gravitat ing test particles that 
are arbitrarily nearby in phase-space ( Valluri & Merritt 2000; 
iKandrup & Sideris 2003). Studies that used live A^-body sim- 
ulations (e.g., Kandrup et al. 1994) have compared pairs of or- 
bits in two different realizations of simulations where the ini- 
tial conditions had been perturbed by an infinitesimal amount. 
Here we work with orbits in a given A^-body simulation and 
are therefore limited (by the resolution of the simulations) in 
our ability to choose pairs of particles that are very close in 
phase-space. 

We select a large number (1000) of nearest-neighbor pairs 
in one of the two merging A^-body halos and follow their sep- 
arations through the entire merger process. We have checked 
that when two identical halos merge, the results were inde- 
pendent of the halo chosen to perform the analysis. To pick 
particle pairs we adopt the following scheme: 

1 . All the particles in a given halo are sorted in their sep- 
aration r from the MBP of that halo, and the particles 
are binned in ten radial shells of equal width, extend- 
ing from the MBP to the initial virial radius of the halo 
(rvii - 256.7 kpc). In the analysis that follows we com- 
pare behavior of particles in the 1st, 4th and 7th shell 
from the center of the potential. The outer radius of the 
first shell is just outside the scale radius, the 4th shell 
lies at the half mass radius and 7th shell lies at the 3/4 
mass radius of the initial NFW halo. 

2. A particle P* with separation rp from the MBP is picked 
at random in a given shell, and the 5000 nearest par- 
ticles in configuration space to the particle P* are se- 
lected. 

3. The velocity separations |v| of each of the 5000 parti- 
cles relative to particle P* are computed, and the 2500 
particles with the smallest |v| separation from P* are se- 
lected. 

4. The magnitudes of the separation of total angular mo- 
mentum J of each of the 2500 particles above relative 
to particle P* are computed, and the 1250 particles with 
the smallest J separation from P* are selected. 

5 . The angle between the angular momentum J and the an- 
gular momentum Jp. is computed for each of the 1250 
particles, and the 625 particles with the smallest angle 
are selected. 

6. From the remaining 625 particles, the particle with the 
total energy E - T + V closest to the particle P* is then 
picked as its nearest neighbor. 



In general, after an initial phase of exponential growth, the 
separation of two nearby particles saturates and then oscil- 
lates about some median value, due to small but finite differ- 
ences in th e initial oscillation frequencies in the global poten- 
tial (e.g.. IKandrup & Siderisll20b3l) . To determine the behav- 
ior of a "typical pair" of particles in a given shell, we first 
pick 1000 pairs of particles in each shell. We then report the 
median value of the separation of that phase-space quantity 
as a function of time for the 1000 particle pairs. We find (as 
is well known) that the median is a more reliable estimator 
than the mean separation of particles because each particle 
has a small but finite probability of being ejected from the 
system as a result of the merger, and such ejections cause a 
small fraction of particles to form a long tail in the histogram 
of separations at each time step. Thus, for all quantities de- 
fined below it should be implicitly understood that we report 
the median value of the parameter for 1000 pairs of particles. 
Specifically, we monitor the following quantities. 

1 . "In(Ar)" - the natural logarithm of the separation in the 
spatial coordinate r of the median pair of particles. 

2. "ln(|Av|)" - the natural logarithm of the absolute value 
of the relative three dimensional velocity v of the me- 
dian pair of particles. 

3. "A/i" - the separation in total energy of the median pair 
of particles. 

4. "A7" - the separation in the total amplitude of the an- 
gular momentum of the median pair of particles. (We 
pick pairs whose vector directions have angular separa- 
tions of less than 90 deg. As the particles separate, we 
track the separation in the amplitude of J but not the 
direction.) 

The quantities defined above differ in one important aspect 
from previous studies in which similar quantities have been 
used: all particles in the pairs are self-gravitating (and not test 
particles) and are picked from an initial distribution function 
that is self-consistent with the spherical NFW potential. This 
implies that the typical initial pair separation In(Ar) of parti- 
cles depends on their location in the potential such that pairs 
in the central cusp will be much closer than the typical initial 
pair separation in the outer regions. We therefore often scale 
the quantities relative to their values at f = Gyr Moreover, 
especially when calculating the velocity separation In(Av), we 
observe the mutual gravitational attraction of nearby particles. 

The quantities AE and AJ provide additional insights be- 
cause in the case of an isolated A^-body halo, both quantities 
are expected to remain constant (to within numerical errors) 
for the entire duration of the simulation (following the ini- 
tial increase in these quantities due to the Miller instability). 
Since E and J are integrals of motion in equilibrium poten- 
tials, any changes in the separation of these quantities in the 
case of the merger reflects the influence of the time-dependent 
potential. Evolution with time of these quantities for various 
simulations are described in § 13.11 

2.3.2. Mixing of ensembles in phase space 

Another way to quantify the mixing rate is by monitoring 
the secular evolution of ensembles of infinitesimally nearby 
test particles due to the presence of chaotic orbits. We con- 
ducted a series of mixing experiments where we monitored 



6 



Valluri et al. 



the spread of an ensemble of 1000 nearby particles as a func- 
tion of time in two phase-space coordinates. 

Such experiments have been carried out previously to 
study the mixing of microscopic ensembles of chaotic 
orbits in static, non -li near p otentials (jM ahon et al.' 
19951; iMerritt & Valluri| IT996t iKandrup et al. 2000: 
Kandrup & Siopis 2003). These authors picked ensem- 
bles of orbits with initial conditions that uniformly sampled 
an infinitesimally small region of phase-space and measured 
the rate at which the ensemble spreads and fills a region of 
configuration space. Ensembles associated with strongly 
chaotic orbits were found to evolve to a near-invariant distri- 
bution within 5-10 orbital crossing times, while ensembles 
associated with weakly chaotic or "sticky" orbits (orbits that 
lie close to resonance in phase-space) evolved on much longer 
time scales and often did not reach an invariant distribution. 
Ensembles associated with regular orbits only spread slightly, 
due to phase mixing that resulted from the small differences 
in their initial integrals of motion. 

Once again the resolution of the simulations determine how 
closely we can sample the phase-space. We pick ensembles 
of the nearest 1000 particles about a given particle (P*) in the 
higher resolution merger simulation only (run HRBp), follow- 
ing a scheme similar to that for picking pairs of nearby parti- 
cles ( § 12.3.11 ) to obtain the 1000 nearest particles to P*. 

Configuration-space coordinates have traditionally been 
used to study chaotic mixing of microscopic ensembles in 
smooth potentials. In such experiments microscopic ensem- 
bles of regular orbits do not mix at all in configuration space 
so all the observed mixing can be attributed to chaos. How- 
ever tests carried out with ensembles selected by the above 
scheme in the A^-body isolated spherical NEW halo showed 
quite significant spreading in configuration-space. The evolu- 
tion of the ensembles in configuration space coordinates fol- 
lowed the following pattern: ensembles spread rapidly in con- 
figuration space on a short timescale (due to the Miller insta- 
bility, see l2.2l i. and then continued to spread more slowly at a 
rate that scaled with local orbital period. At a given radius the 
mixing rate was found to be sensitive to the particular ensem- 
ble in question - the nature of the parent orbit (whether it is 
box-like or tube-like), the location of the orbit relative to the 
orientation of the merger etc. Since their mixing in configu- 
ration space is dominated by the Miller instability and phase 
mixing, our ensembles are too large to be usefully compared 
with the microscopic ensembles used to measure chaotic mix- 
ing in the works of previous authors. However, we found that 
ensembles selected by the above scheme in the isolated halo 
showed negligible mixing in the coordinates E, J (in agree- 
ment with the previous findings of lHut & Heggiel2001h . Con- 
sequently we focus our attention on mixing in these two vari- 
ables. The main objective of carrying out these experiments 
is to observe mixing in E, J that results from the merger. 

We define the quantities (£, J') for each particle in the en- 
semble as follows: 

E - E(P*) 

^=^(iir 

J-'-^^ (4) 

J(P*) ^ ' 

where E{P*), J{P*) are the energy, magnitude of the total an- 
gular momentum of the particle P* respectively. The quanti- 
ties &,J', are computed for each of the 1000 particles in the 
ensemble. 

The results of these experiments are described in § 13.21 



3. RESULTS 

3.1. Separation of Nearby Particles in Phase-Space 

We begin by demonstrating the ability of In(Ar) to identify 
chaotic behavior. This quantity (and all the others that we use) 
are compared to equivalent quantities measured in an isolated 
equilibrium NEW halo. The three panels of Eigure|2]show the 
initial separation in spatial coordinate In(Ar) (relative to its 
value at f = Gyr) in three different radial shells (1st, 4th and 
7th). The horizontal axis gives time in units of the median 
crossing time t^ in each shell of the corresponding isolated 
NEW halo. The local cross time are = 0.53,2.59,4.95 Gyr 
for shells 1, 4, 7 respectively. We see that in both the merger 
(solid curve) and the isolated halo (dot-dash curve) the initial 
separation is nearly linear in the quantity [In(Ar)/ ln(Ar(0))], 
indicating an exponential change in radial separation. The 
linear increase in In(Ar) on this timescale constitutes a clear 
manifestation of the Miller instability (see § 12.2b . This phase 
of exponential divergence is short lived (< 35% of the lo- 
cal crossing time). The e-folding time of the Miller insta- 
bility t, = 0.39, 1.06,2.73 Gyr for shells 1, 4, 7 respectively 
(roughly half the local crossing time in each shell). In both 
the merger simulation and in the isolated halo, the initial rapid 
change in In(Ar) saturates at the same value - a value deter- 
mined by the length scale over w hich the smooth pote ntial 
dominates over the graininess (Kan drup & S"mMilll991bh . At 
smaller radii, the higher particle densities result in saturation 
values of In(Ar) that are higher that at larger radii. Eigure |2] 
is clear evidence that micro-chaos is detectable despite our 
inability (imposed by our use of a A^-body distribution func- 
tion) to choose particles arbitrarily close to each other in phase 
space. 

We now focus on the behavior of the four quantities (In(Ar), 
ln(A|v|), IS.E and A/) over longer time-scales. As before all 
quantities are compared with the equivalent quantities mea- 
sured in an isolated equilibrium NEW halo to control for ex- 
ponential Miller instability. Eigure [3] shows the evolution 
of In(Ar) and In(Ay) as a function of time in 3 different ra- 
dial shells (1st, 4th and 7th from the center). The quantities 
plotted are scaled using the initial values of the quantities at 
t - Q Gyr; i.e., the top panels show [In(Ar)/ ln(Ar(0))] and the 
lower panels show [ln(A|v|)/ ln(A|v(0)|)] (where Ar(0),Av(0) 
are the separations at f - 0). The solid lines show the evolu- 
tion of the quantities in the low-resolution merger simulations 
(run Bp 1 : two NEW halos with 2x10^ particles per halo), and 
the dot-dashed curves show the evolution of these quantities 
in the isolated equilibrium A^-body NEW halo. In what fol- 
lows we point out some characteristic features of this figure 
that are seen elsewhere in this section. 

The lower panels of Eigure[3]show the separation in ln(A| v|) 
in both the isolated halo and the merging halos. There is a 
noticeable initial decrease in ln(A|v|) (which occurs simulta- 
neously with the linear increase in [In(Ar)/ ln(Ar(0))]), which 
we attribute to the mutual gravitational attraction between the 
pair of particles that causes a deceleration of the particles 
while their radial separation increases. The initial decrease 
in ln(A|v|) is seen in both the merging halos and the isolated 
halo, confirming that this is not due to the merger but is a 
characteristic of the A^-body simulation. 

In the innermost (1st) shell, after the initial growth in In(Ar) 
and the initial decrease in ln(A|v|), there are periods of time 
during which these quantities are essentially constant, sep- 
arated by several sharp decreases in In(Ar) and correspond- 
ingly sharp peaks in ln(A|v|) not seen in the isolated halo. The 
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Fig. 2. — Divergence in the quantities In(Ar)/ ln(Ar(( = 0)) as a function of time (in units of local crossing time) for particles in the 1st, 4th and 7th radial shells. 
The solid lines are for the merger of two NFW halos (run Bpl). The dot-dashed lines are for the control simulation of an isolated NFW halo. The dashed straight 
lines are best-fits to both the solid lines and the dot-dashed lines in the region tjtc = [0.05 - 0.35]. 
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Fig. 3. — The top three panels show divergence in the quantities In(Ar)/ ln(Ar(f = 0)) as a function of time for three different shells in the halo (from left to 
right: the first, fourth, and seventh shells respectively). The lower three panels show the quantity ln(A|i'|)/ ln(A|v(? = 0)|), for the same three shells. The solid 
lines are for the merger of two NFW halos (run Bpl). The dot-dashed lines are for the control simulation of an isolated NFW halo. 



dips in In(Ar) and peaks in ln(A|v|) are coincident with each 
other and correspond to the pericenter passages of the MBPs 
of the two halos seen in Figure [1] It is the first and second 
pericenter passages that cause the most significant changes in 
the particle separations (In(Ar), ln(A|v|)). 

It is quite remarkable how little systematic change in sepa- 
ration of either r or v is seen in between pericenter passages, 
suggesting that very little change in the macroscopic separa- 
tions of particles occurring at any time during the merger ex- 
cept during the pericenter passages. This can be understood as 
follows: during pericenter passages of the centers of masses 
of the two halos, there is significant overlap in the two po- 
tentials leading to a sudden deepening in the potential experi- 
enced by the particles leading to a transient compressive tidal 
field that simultaneously decreases the separations in parti- 
cles (seen by the dips in In(Ar)) and increases their kinetic 
energies (seen in the sharp increases in ln(A|v|)). Dynamical 



friction between the two halos is also strongest at pericenter 
and consequently the majority of their orbital energy and an- 
gular momentum are lost impulsively at these times. When 
the two halos separate enough that their innermost shells no 
longer overlap they experience little or no change in their po- 
tentials, so there is little further change in separation of tra- 
jectories. Thus, the global potential fluctuations that occur 
between pericenter passages do not cause any further increase 
in particle separations. Behavior similar to that in shell 1 is 
seen in shells 2 and 3 but is not shown here. 

The effect of the first pericenter passage is seen to occur 
simultaneously at all radii. In shells 4, 7 we observe a dip 
in In(Ar) and peak in ln(A|v|) at the first pericenter passage 
which is both broader and less intense than in shell 1. At 
larger radii the first pericenter passage is followed by a con- 
tinuous change in In(Ar) and ln(A|v|) but subsequent pericen- 
ter passages are less clearly visible. This indicates that while 
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there is no propagation delay in the impulsive tidal shock with 
radius, the response is longer lived at larger radii. This is be- 
cause, the centers of the two merging halos do not separate 
beyond the radius of the 4th and 7th shells after the first peri- 
center. Consequently the change in the external potential is 
not as sudden as for the inner shells. 

As discussed in § l2.3l it is extremely difficult to make a re- 
liable quantitative measurement of an exponential instability 
on a short timescale. However, both the qualitative and quan- 
titative behavior of initially nearby particles seen in Figure 2 
indicates that we are in fact detecting the Miller instability. If 
the oscillating external potential is also contributing to an ex- 
ponential instability of orbits one might expect that if pairs of 
nearby particles were picked during the phase when the po- 
tential is changing very rapidly, the rate of separation in some 
quantity (say r) would be greater than in either the isolated 
halo or at the start of the merger Figure|4]compares the sepa- 
rations of particles In(Ar) when pairs were picked at f = Gyr 
with pairs picked at f = Gyr in the isolated halo and sepa- 
ration of pairs of particles picked at f = 3 Gyr. For the pairs 
of particles picked at f = 3 Gyr in the merger, the time axis 
is f - 3 so that initial slopes of the three curves can be more 
readily compared. We chose a start time t - 3 Gyr to compare 
against since it corresponds to a time at which the virial ratio 
is changing very rapidly and corresponds to the first apocenter 
separation of the MBPs of the two halos (identical behavior 
was found when particles were picked at f = 4 Gyr). 

Since the potential is changing rapidly at this time, we 
might expect to observe an enhanced rate of separation of 
nearby pairs of orbits due to chaos induced by the potential 
fluctuations provided the e-folding time for this instability is 
comparable to or greater than that due to the Miller instability. 
This figure shows, on the contrary, that the initial exponential 
separation of particles in r at f = 3 Gyr does not have a sys- 
tematically greater slope than the rate of separation of parti- 
cles at f = Gyr In shell 1 the slope of the dot-dash curve 
is smaller than that of the solid curve (which we attribute to 
the fact that heating from the first pericenter passage caused 
nearby particles to be further apart in phase space); in shell 
4 the slopes of all three curves are essentially identical; only 
in shell 7 is the slope of the dot-dash curve slightly steeper 
than that of the black curves indicating that chaoticity may be 
playing a bigger role at large radii. 

While this is clearly not conclusive evidence, it is sug- 
gestive of the hypothesis that chaoticity due to the time- 
dependent potential, is not significant compared to the 
chaoticity due to the Miller instability. It may play a small 
role in the evolution of the system at larger radii. However 
there is little evidence to suggest that chaotic mixing from or- 
bits in the time-dependent potential is driving the relaxation. 

To better understand the mechanism driving the relaxation 
we now look at separation of particles in energy and angular 
momentum space. In integrable time-independent potentials, 
energy and some quantity akin to angular momentum are gen- 
erally integrals of motion. It is therefore of interest to quantify 
how these quantities change during a merger. 

Figure|5]shows the separations in energies (top panels) and 
angular momentum^ (bottom panels) for the same pairs of 
particles as in Figure |3] As before, the dot-dashed lines are 
for pairs of particles in the isolated halo, while the solid lines 
are for pairs of particles in one of the merging halos. The 

' The angular momentum (J) is about an axis perpendicular to the plane 
containing the orbits of the center-of-masses of the two merging halos. 



quantities AE, AJ are scaled relative to their initial values at 
t - Gyr. AE and AJ (as expected) are very close to zero in 
the isolated halo** .The sharp initial rise in AE and AJ reflects 
the response of the particles to the presence of second halo. 
This a consequence of the fact that although each individual 
halo is in equilibrium separately, when the two halos are set 
up on a parabolic trajectory at f = they are no longer in equi- 
librium, since particles now experience the external potential 
of the second halo. 

After the initial increase in AE and AJ, there are sharp but 
transient increases in the separation of these quantities that 
also occurs primarily during pericenter passages. In fact, the 
peaks in AE and AJ are strongly correlated with the pericen- 
ter passages in all radial shells and provide striking evidence 
that it is these events that are primarily responsible for scat- 
tering particles in E and J. We also note that the step-like 
changes in these quantities and their correlation with pericen- 
ter passages is seen at all radii unlike in the case of Ar and 
Av which changed more slowly and continuously in the outer 
regions. In the inner most shells there is little or no sepa- 
ration in phase-space quantities occurring between pericenter 
passages. The slight increase in AE and AJ seen in shell 4 
and shell 7 between f = 2 Gyr (first pericenter passage) and 
f = 5 Gyr (second pericenter passage) could be evidence for 
slower chaotic mixing. After the first pericenter passage the 
outer regions of the two halos never completely separate again 
and the continued overlap of the two potentials means that 
at larger radii the tidal compressive field and dynamical fric- 
tion is able to operate over a longer duration and consequently 
more gradual changes occur as well. It is clear, however that 
the impulsive changes in AE and AJ during multiple pericen- 
ter passages have the strongest effect on the mixing in E and 
J 

Figure|6]compares separations in all four quantities In(Ar), 
ln(A|v|), AE and AJ for two merger simulations with different 
numerical resolutions. We recall that the run HRBp was per- 
formed with a factor of 10 more particles and approximately 
a factor of 2 higher force resolution than the run Bpl. We 
present the comparison for only the innermost shell since this 
is expected to be the most susceptible to numerical effects, if 
any. There is clearly very little difference between the overall 
behavior of separations of particles in any of the 4 quantities 
plotted. The only noticeable difference is in the initial behav- 
ior of In(Ar), ln(A|v|). 

The insert in the top most panel shows that the initial radial 
separation of particles in run HRBp is smaller at f = Gyr 
than it is for run Bpl . This is purely a consequence of the fact 
that by increasing the mass resolution the phase-space coor- 
dinates become more densely sampled and this allows us to 
select pairs of particles that are closer in phase-space (see in- 
serts in top two panels). The insert in the second panel shows 
that the initial velocity separation (ln(A|v|)) decreases to lower 
values in run HRBp compared to run Bpl. This also points to 
the increased gravitational deceleration due to the proximity 
of the nearest neighbor and the decrease in the softening pa- 
rameter. 

All curves are plotted without scaling them relative to their 
initial separations to illustrate that apart from differences in 
the initial separation (which we attribute to the Miller insta- 

We note that they are not precisely zero because, although the total en- 
ergy of the system is conserved to within numerical error, the computation of 
(AE/AE(t = 0)), E for a given particle is relative to the MBP - which can be 
a different particle at each time-step. 
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Fig. 4. — The panels show divergence in the quantities ln(Ar)/ ln(Ar(? = 0)) as a function of time for 1st. 4th, 7th shells. The solid curves are for the merger of 
two NFW halos (run Bpl) with pairs of particles selected axt = 0. Gyr. The dashed curves are for pairs of particles selected at f = 3.0 Gyr in the isolated halo and 
the dot-dashed curves are for pairs of particles selected at t = 3.0Gyr. The curves for the pairs starting aX t = 3.0 Gyr have been translated in time to ? = Gyr 
The dot-dashed curves and the dashed curves indicate that the initial exponential instability of nearby orbits during the merger is not very different from that in 
the isolated halo, independent of the time at which the particle pairs are selected. 
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Fig. 5. — The top three panels show divergence in the quantities AE/AE{t = 0) as a function of time for three different shells in the halo (from left to right: the 
first, fourth and seventh shells respectively). The lower three panels show the quantity AJ/AJ(f = 0), for the same three shells. The sohd lines are for the merger 
of two NFW halos (run Bpl). The dot-dashed lines are for the isolated NFW halo. 



bility), the separations in all four quantities in the different 
resolution mergers saturate at the same values. This is con- 
firmation that particle separations saturate at values that are 
determined by the global dynamics of the merger and not by 
the resolution of the simulations. 

Figure |7] shows separations in In(Ar) and In(Av) for parti- 
cles in different radial shells in run Bpl and in the merger of 
an NFW halo with a shallow cusp. The behavior is remark- 
ably similar in both cases except in the innermost shell. This 
is expected, since the innermost shell is the only one in which 
the density profiles (and consequently the phase-space density 
distributions) differ significantly. It is well known from pre- 



vious studies dBovlan-Kolchin & Mall2004l : iKazantzidis et al.l 
IJOOe") that cusps in shallow potentials are less robust than 
cusps in NFW halos. Particles in the inner most shell of the 
shallow cusp experience a larger increase in energy due to 
each pericenter passage than particles in the NFW halo be- 
cause they experience a greater relative increase in the depth 
of the external potential arising from the overlap of the two 
halos, and consequently the phase space volume accessible to 
them is larger Beyond the scale radius the two halos have es- 
sentially identical density profiles, and separations in In(Ar) 
and In(Av) are almost indistinguishable. 
The results of our investigation of particle separations in 
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Fig. 6. — The four panels (top to bottom) show divergence in the four 
quantities In(Ar), In(A|v|), A£ and AJ plotted as a function of time for pairs 
of particles in the innermost shell of one of the two merging NFW halos. 
The solid lines are for the merger of two halos with 2 X lO' particles in each 
halo (run Bpl); dot-dashed lines are for the merger of two halos with 2 X 10^ 
particles in each halo (run HRBp). The inserts in the top two panels show 
differences in the initial separation due to the Miller instability. The dot- 
dash line starts at lower values in both In(Ar), ln(A|i'|) which is a result of 
picking particles closer in phase space in ran HRBp. For ran HRBp \n(Ar) 
has a shorter e-folding time, due to the higher particle number and smaller 
softening length. Apart from the behavior of initial separations in all four 
quantities, the low and high resolution mergers are almost identical. 

In(Ar), ln(A|v|), AE and AJ can be summarized as follows: 

1 . The initial separation of pairs of nearby particles in A^- 
body systems such as those studied here is a result of 
the exponential instability of the A^-body problem, the 
so called "Miller instability." The qualitative and quan- 
titative behavior of the separation of particles: e.g. their 
dependence on local crossing time (Fig. O, their de- 
pendence on softening and force resolution (Fig.|6]l are 
completely consistent with previous studies of this in- 
stability. It is important to point out that the very fact 
that we are able to detect the exponential divergence of 
orbits due to the Miller instability indicates that despite 
the fact that pairs of nearby particles in the A^-body sim- 
ulations are separated by macroscopic distances (i.e., 
are not infinitesimally close by), our method for identi- 
fying an exponential instability can be applied success- 
fully to an A^-body system. 



2. Subsequent to the saturation of the separations due to 
the Miller instability, the most significant increase in 
separations of nearby particles (in r, v, E, or J) occur 
during pericenter passages of the MBP of the two halos 
and is a consequence of the compressive tidal shock- 
ing (due to the overlap of the two halos) and dynam- 
ical friction between the two halos that occurs during 
the pericenter passages. In the innermost shells there 
is little or no separation in phase-space quantities oc- 
curring between pericenter passages, despite the large 
global fluctuations in potential occurring during these 
times. At larger radii separation in r, v continue follow- 
ing pericenter passages. At larger radii, there is also a 
small (but short lived) increase in separation in E and J 
between the first and second pericenter passages. 

3. If we consider the increase in the macroscopic separa- 
tion of particles in (r, v) to be a measure of the mix- 
ing of particles in these quantities, these plots would 
lead us to conclude that the majority of mixing occurs 
due to the Miller instability and several phases of mix- 
ing that occur during pericenter passages of the two ha- 
los. In all radial shells of merger simulation we observe 
that the macroscopic separation of initially nearby par- 
ticles increases until it saturates. In the merging halos, 
saturation occurs at larger values of In(Ar), ln(A|v|) re- 
spectively, than in the isolated halo. This is a conse- 
quence of the fact that the accessible phase-space vol- 
ume has increased during the merger process. We saw 
from Figure |5] that AE and AJ (which represent in- 
creases in the available phase-space states accessible to 
particles) increased in step-wise fashion primarily dur- 
in g pericenter passages of the MPBs, as was predicted 
bv lSpergel & Hernquisti (119921) . By f = 7 Gyr the virial 
ratio 2r/|y| 1, although the potential continues to ex- 
perience long-lived but low-amplitude oscillations until 
13 Gyr. The low-level potential fluctuations cause little 
change in the separation of either r or v beyond 7 Gyr. 

In Lynden-Bell's (1967) model for "violent relaxation", 
particles experience changes in their phase-space coordinates 
due to their interaction with a time-varying background po- 
tential. It has been argued by various authors that a time- 
dependent potential alone is inadequate to cause relaxation, 
since nearby particles in phase space will merely be rela- 
beled in energy and will not separate (or mix). The results 
from this section indicate that during the merger of two ha- 
los, nearest neighbor particles in phase space do mix in E 
and J . Further, the primary causes of mixing are the com- 
pressive tidal shocks and dynamical friction which transfer 
energy and angular momentum from the orbital motion of the 
two m erging halos to the p articles during the pericenter pas- 
sages (ISpitzer & Chevalie r 1973; Gnedin et al. 1999; Vallur^ 
Il993h . The particles respond to sudden increments in their 
orbital energy and angular momentum by jumping to new or- 
bits, resulting in mixing in phase space and evolution to a new 
distribution. This latter phase of evolution could indeed be the 
result of some weaker form of chaotic mixing, coupled with 
phase mixing, both of which occur on a longer timescale. 

3.2. Mixing of ensembles in phase-space 

In this section we present results of the mixing experi- 
ments with an initially coherent ensemble of particles selected 
from the A^-body simulation. As discussed in 12. 3. 21 we con- 
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Fig. 7. — Solid lines show separations in In(Ar) and In(Av) for the merger of two NFW halos and the dot-dash lines show the separations of these quantities for 
the merger between an NFW halo and a halo with an inner shallow cusp of y = 0.2. The behavior is similar in both cases except in the innermost shell where the 
differences in the density profiles are significant. 



fine our analysis to mixing in the scaled energy and angular- 
momentum variables £, ST. 

Figure [8] shows the evolution of an ensemble of 1000 self- 
gravitating particles in the 4th shell in run HRBp. The density 
contours are obtained using a kernel smoothing algorithm. At 
f = 2 and f = 5 Gyr which correspond to the first and sec- 
ond pericenter passage, respectively, the ensemble undergoes 
a sudden spreading in energy (£) which is a consequence of 
the deepening of the potential during the overlap of the two 
halos. By f = 6 Gyr the particle distribution fills a triangular 
region in £, ST space. A greater spreading in £ at the small- 
est values of angular momentum is a consequence of the fact 
that these orbits are on the most radial orbits and consequently 
experience greater changes in their potential energies. Parti- 
cles with the smallest (most negative) values of £ experience 
the greatest spread in angular momentum largely as a conse- 
quence of their being ejected during the tidal shocks to larger 
radii. 

We carried out similar mixing experiments for over 25 dif- 
ferent ensembles in each shell. The behavior of the ensemble 
in Figure[8]was found to be quite representative of all the dif- 
ferent ensembles. 

Figure |9] shows contours of projected density of particles 
as a function of the quantities {&,J~) in run HRBp, for en- 
sembles in 3 different radial shells at 4 different times in the 
evolution. The evolution of the ensembles in shells 1, 4 and 7 
are very similar. Ait = 2 Gyr (first pericenter passage) the en- 
semble spreads in £. Similar increases in spread and are seen 
during the 2nd pericenter passages (but are not shown here). 
After the first pericenter passage the majority of the particles 
in the ensemble at t - 4 Gyr return to a more compact dis- 
tribution in (£, but one with a larger spread in ST at small 
£. The bottom panels correspond to f = 6 Gyr, after the first 
three pericenter passages have occurred. At this time all the 
ensembles fill a roughly triangular region in £, J' space. This 



characteristic triangular final distribution is seen in all shells 
and for essentially all of over a hundred different ensembles 
we examined. This final distribution is not uniform in density, 
but the density contrast across the ensemble is much smaller 
than in the initial distribution. We note that when these en- 
sembles were re-observed at 15 Gyr there was slightly greater 
uniformity of density in the lower-density tails but otherwise 
they had changed very little from their distributions at 8 Gyr 

There are two important observations that can be made: (a) 
changes in £, J' occur primarily during pericenter passages 
(Fig. [8]l , and (b) the final distribution in S,J' space of the 
ensembles at all radii are similar, and the rate at which this 
final distribution is reached is independent of radius (Fig.|9]). 

The median crossing time in the innermost shell (r < 
25 kpc) is ~ 5 X 10** yrs, while in the 7th shell is ~ 5 Gyr 
It is obvious from Figure|9]that although the crossing time in 
shell 7 is a factor of 10 longer than in shell 1, the diffusion 
of the ensembles in (£, ^-space appears to occur primarily 
due to the pericenter passages of the MBPs and there is only a 
slightly increased rate of spreading at larger radii. Phase mix- 
ing and chaotic mixing in isolated potentials occur at rates that 
scale with the local dynamical time so one would expect that 
mixing effects that depend on the local dynamical time (such 
as chaotic mixing) would occur faster at small radii than at 
large radii. The fact that the evolution of the ensembles at dif- 
ferent radii occur episodically following pericenter passages, 
rather than continuously at a rate that scales with local orbital 
time points to the importance of the impulsive tidal events as 
the processes that drives the transition to an equilibrium dis- 
tribution in energy and angular momentum. 

Indeed, the absence of a strong dependence on local cross- 
ing time in Figure|9]provides the strongest evidence so far that 
mixing in energy and angular momentum is driven by com- 
pressive shocking and dynamical friction that occur at peri- 
center passages. While chaotic mixing, as defined for static 
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Fig. 8. — The evolution of an ensemble of 1000 nearest neighbors in phase-space in shell 4 at < = (run HRBp) plotted at 1 Gyr time intervals in the merger The 
plot shows isodensity contours for the particle distribution. The red contours correspond to the highest density regions and black to the lowest density regions. 
The contours are spaced at logarithmic density intervals relative to the maximum density contour to enhance the visibihty of low density regions. Pericenter 
passages (at ^ = 2 and t = 5 Gyr) are seen to cause a sudden spreading of the entire ensemble in £ (a result of the tidal shocking at these times). Two additional 
pericenter passages occur between t = Sandt = 6 (see Fig. 1). By f = 6 Gyr the ensemble has experienced 4 pericenter passages and has settled to fill a triangular 
region in the £, J" space. 



smooth potentials, might be occurring during the merger of 
two halos, it plays a minor role in driving the remnant to equi- 
librium. 

4. INCOMPLETE RELAXATION AND 
REDISTRIBUTION OF PARTICLES IN MERGERS 

Two interesting aspects of merger simulations of DM ha- 
los that have been pointed out recently are (a) that cusps 
are remarkably robust and (b) that about 40% of the to- 
tal mass of the remnan t lies beyond fiducial virial radius 
jKazantzidis et al.ll2006l) . Two questions that arise from these 
findings are: 1) how do the final radial distributions, ener- 
gies and angular momenta of particles in the remnant depend 
on their original location in the merging halos and, 2) From 
where do the particles that end up outside the virial radius 
in the merger remnant originate? Addressing the first ques- 
tion is important since the merger causes the redistribution of 
particles to occur in a way that the density profiles and phase- 
space distribution functions preserve homology. This may be 
interpreted as an indication that particles at all radii are being 
heated uniformly at all radii. However, it is also reasonable to 



assume that less bound particles in the outer part of the halo 
are preferentially heated. 

It is also of interest to understand how the relative change in 
energy or angular momentum of particles depends on their ini- 
tial location in the halo. In this section we examine the redis- 
tribution of particles in radius, energy and angular momentum 
and the dependence of their final locations on their initial lo- 
cation in the merging halos. We note in passing that while the 
isolated halos and the initial merging halos are spherical, jus- 
tifying the use of spherical radial shells, the merger remnant 
is significantly prolate-triaxial with minor- to-major principal 
axis ratio c/a varying from 0.5 at the center to 0.7 at the virial 
radius. Although the merger remnants exhibit significant de- 
partures from spherical symmetry, in what follows we ignore 
the triaxiality of the final distribution when binning particles 
in radius. 

Each panel in Figure[TO]shows how particles that were orig- 
inally in a given shell (indicated by the caption) are redis- 
tributed in the final radial shells at f = 9 Gyr when the merger 
would conventionally be deemed complete. The results for 
isolated halo are an important baseline for comparison. 
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Fig. 9. — Isodensity contours for ensembles of particles plotted as a function of £, J for shell 1 (left hand column), shell 4(middle column) and shell 7(right-hand 
column) plotted for dilferent times in their evolution (see text for details). 



In the top panel (shell 1), the solid lines (merger) and dot- 
dash lines (isolated halo) are remarkably similar but indicate 
that while the cusp is quite robust, in fact as much as 20% of 
the particles in the isolated halo and over 40% of the particles 
in the merging halos end up outside the cusp at the end of the 
simulation. The primary reason is that 80% of the particles in 
the isolated halo have apocenters within the cusp. The merger 
only causes a small additional fraction (20%) of these bound 
particles to be spread to shells immediately outside the cusp 
(shells 2-4). In comparison with the outer shells, this level of 
redistribution for shell 1 is somewhat modest. 

The middle panel (shell 4) shows that, for the isolated NFW 



halo (dot-dashed curve), ~ 50% of the particles are in shells 
3 and 4 at the end of the simulation, but the remainder are 
almost uniformly redistributed within shells 2, and 5-10. In 
shell 7 of the isolated halo ~ 45% of the particles stay in 
shells 6 and 7 while the remainder are redistributed almost 
uniformly to all radii from shell 4 outward. Thus, there is 
significant radial redistribution of particles even in the non- 
evolving isolated NFW halo purely because the initial dis- 
tribution function has an isotropic velocity distribution with 
roughly equal fractions of radial and tangential orbits. The 
peaks in shells (3,4) and (6,7) are due to the more tangen- 
tially biased orbits. A comparison of the dot-dash curves with 
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Fig. 10. — The fraction of particles (W/Wmax) that lie in each of the radial 
shells at the final time step. The vertical dashed line and the caption in each 
panel indicate the shell in which particles were at ( = Gyr. Particles from 
one of the two NFW halos (in run Bpl) is shown as a solid line while particles 
in the same shells of the isolated NFW halo are shown by dot-dashed Une. 



the solid curves shows that the main effect of the merger is 
to flatten out the peaks due to tangential orbits thus making 
the final distributions more radially anisotropic. At all radii 
the merger is responsible for redistributing about ~ 15 - 20% 
more particles from a given bin to other radial bins. 

A similar study of the merger of a NFW halo with a halo 
having a shallow cusp showed that except in the innermost 
shell where the profiles of the two halos differ, the particles 
in the two halos with cusps of different slope are redistributed 
identically. A larger fraction (~ 45%) of particles in the shal- 
low cusp were redistributed by heating to radial shells (2-4) 
directly outside the cusp (in comparison with 20% in NFW- 
NFW halo mergers). The differences in the distribution of 
particles at the end of the merger confirm that NFW cusps are 
robust and retain a much more significant fraction of their par- 
ticles than do shallow cusps, c onfirming what has been pre- 
viously found by other authors dBovlan-Kolchin & Mall2004l : 
iKazantzidi s et alj|2006l) . 

Figure [TT| shows distributions of the changes in the kinetic, 
potential, and total energies of all particles in different shells 
at the end of the merger (A/i = £finai-£^initiai)- The y-axis gives 
the fractional change in number of particles (relative to the to- 
tal number of particles in a given radial shell). Particles in the 
innermost shell experience a net decrease in total energy that 
is a result of the overall decrease in the potential energy of the 
particles poi nting to an energy segregation phenomenon dur- 
ing mergers dFunato et aljl992h . This is expected, since there 
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Fig. 1 1 . — Histograms of the change in energy (Ail = -Efi„al ~ ^^iniiial) for all 
particles in shells 1 (top), 4 (middle), and 7 (bottom). The solid line shows 
the change in total energy, the dot-dashed shows the change in kinetic energy, 
and the dashed line shows the change in potential energy. 



is an overall increase of ~ 60% in the mass of the cusp com- 
pared to the mass in the initial cusp (Kazantzidis et al. 2006). 
At all other radii, there is a net increase in total energy that 
is most significant at larger radii. This is a result of signifi- 
cant fractions of particles gaining potential energy and being 
redistributed to larger radii as seen from the multiple peaks in 
the potential energy distributions arising from particles heated 
during each pericenter passage. It is striking to note that in 
all shells the kinetic energy distributions remain peaked about 
A/i = and are more peaked than Gaussian, indicating that 
the majority of particles experience only a small change in 
their kinetic energies despite experiencing large changes in 
their potential energies. This is additional confirmation that 
the redistribution in energy that we saw in the mixing experi- 
ments 03.2b is primarily due to changes in potential energies 
of particles due to interactions with the background potential 
during pericenter passages. 

The innermost shell includes all particles within ~ 25 kpc, 
which is slightly larger than the initial scale radius (r^ - 
21 kpc) of the merging NFW halos. Figure [12] is similar to 
Figure [TT|but shows how particles in three equal radial shells 
within r, are redistributed in energy. It is clear that the parti- 
cles in the inner 33% of the scale radius experience the great- 
est deepening in the potential. 

Figure[T3]shows histograms of change in (the component 
of J perpendicular to the orbital plane of the merger) for all 
particles in radial shells 1, 4 and 7. The highly peaked distri- 
butions in the innermost shell indicates that there is virtually 
no net change in for the particles here. In the outer two 
shells, M z has an increasingly wider distribution. Although 
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Fig. 12. — Histograms of the change in energy {AE = iTunai - ^initial) for all 
particles in three equal radius bins within the scale radius (r, = 21 kpc). The 
solid line represents total energy, the dashed-dot represents kinetic energy 
and the dashed represents potential energy. 



it is not plotted here, we note that in the outermost shells 
(shells 9, 10) the median of the distribution of AJ^ shifts to 
slightly positive values, indicating a small gain in net angular 
momentum for the outermost particles. The observed increase 
is a result of transfer of the orbital angular momentum of the 
merging halos into the angular momentum of individual par- 
ticles. Our results indicate that the angular momentum of the 
merger is absorbed primarily by particles at large radii. 

As was noted previously, about 10% of the mass of the 
initial DM halos lies outside the fiducial virial radius, but 
it has been found that nearly 40% of the mass of the final 
remnant lies outside the virial radius of the remnant halo 
feazantzidis et al.ll2006l) . We now examine the distribution 
of particles in the merger remnant of run Bpl to determine 
where these particles originate from. Figure [14] shows the 
fraction (Fyir) of the total number of particles lying beyond 
the virial radius of the remnant that originated from each of 
the 10 shells of the original halos. All particles that lay out- 
side the virial radius in the original halos are assigned to shell 
1 1 (which extends from the virial radius to the outer edge of 
the simulation volume). The highest fraction (about 25%) of 
the particles outside the virial radius of the remnant were al- 
ready outside the virial radius of the initial systems. The in- 
nermost shell is the most robust with fewer than (< 1%) be- 
ing ejected beyond /?vir. Interestingly, all shells from the half 
mass radius onward (shell 3 and beyond) contribute roughly 
equally (~ 8 - 1 1%) to the particles that lie outside /?vir in the 
final remnant. This figure and Figure [TO] together show that 
the radial redistribution of particles occurs essentially inde- 
pendently of the original radius of the particle, provided that 
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Fig. 13. — Histograms of change in (component of J perpendicular to the 
orbital plane of the merger) for all particles in the selected shells. 



the particles lie outside the central 3 shells. 

These results further highlight the fact that mass, in its stan- 
dard virial definition, is not additive in mergers. Models in 
which the mergers double the mass within the virial radius of 
the remnant greatly overestimate both the mass and density of 
the merger product. 

5. SUMMARY AND DISCUSSION 

The principal goal of this study was to investigate the 
processes that drive the evolution of self-gravitating systems 
to equilibrium and cause mixing in phase-space. The goal 
was also to determine whether the possible presence of large 
fractions of chaotic orbits arising from the time-dependent 
potential was responsible for driving the system to equilib- 
rium. We summarize the main results below. 



1. All orbits in the A^-body simulation exhibit the well- 
documented "Miller instability" of the A^-body prob- 
lem. This instability results in an initial nearly exponen- 
tial separation of nearby orbits in configuration space 
and energy, and a decrease in velocity separation. The 
e-folding time (f^,) for the separation increases with in- 
creasing softening length and decreasing particle num- 
ber. The initial exponential separation saturates on a 
short timescale (~ 0.4fc)- The e-folding time and the 
number of e-folds after which the instal)ility saturates 
are identical in the merging and isolated halos indicat- 
ing that this is purely a consequence of the A^-body na- 
ture of the problem and not a result of the merger. To 
this end, we use control simulations of isolated halos as 
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Fig. 14. — Fraction of particles that originated in each of the radial shells 
(1-11) in one of the original NFW halos involved in the merger that end up 
outside the fiducial virial radius of the remnant. (All particles beyond the 
virial radius are assigned to shell 11.) 

a baseline for comparison with merging halos to iden- 
tify the mixing processes occurring in mergers. 

2. After the initial growth in In(Ar) and decrease in 
ln(A|v|), there are periods of time for which these quan- 
tities are almost constant, separated by several sharp 
decreases in In(Ar) and correspondingly sharp peaks in 
ln(A|v|) that are not seen in the isolated halo. The tran- 
sient changes in In(Ar) and ln(A|v|) are coincident with 
each other and correspond to the pericenter passages of 
the centers of the merging halos. Pericenter passages 
also cause sharp transient step-like changes in the en- 
ergy and angular momentum separation of pairs of par- 
ticles. Between pericentric passages there is very little 
evidence for for mixing in the innermost shell, but at 
larger radiiln(Ar) and ln(A|v|) continue to grow slowly 
with time due to the overlap of the merging halos. 

3. The changes in AE and AJ are more strongly correlated 
with pericenter passages at all radii pointing to the im- 
portance of compressive tidal shocks. At larger radii 
there are is also some growth in separation between 
pericentric passages, possibly pointing to the contin- 
ued action of chaotic mixing. If we assume that sep- 
arations of nearby particles in phase-space are linked 
to "mixing", this implies that the principal drivers of 
strong mixing are the compressive tidal shock and in- 
creased dynamical friction that occurs during the peri- 
centric passages. 

4. Mixing experiments with ensembles of 1000 nearest 
particles in phase-space showed that mixing in the 
phase-space variables £, J' also occurs primarily during 
pericenter passages. Mixing leads ensembles of nearby 
self-gravitating particles to spread and fill a triangular 
region in £, J' space. Mixing in these parameters oc- 
curs at the same rate at all radii, largely independent of 
the characteristic crossing time of the ensemble. This 
is further evidence that the mixing process is driven 
by pericentric passages and is unrelated to mixing pro- 
cesses that scale with the local crossing time. 

5. We confirm the findings of several other authors that 
cusps of DM halos are remarkably robust. The robust- 
ness of cusps is a consequence of the fact that 80% of 
the particles in the first shell have apocenters that lie 
within the first shell. During the merger only a small 
fraction (20%) of particles with apocenters within the 



original cusp are ejected to larger radii. The majority 
of the ejected particles do not get beyond shell 4. A 
much larger fraction (~ 45%) of particles in the central 
regions of core-like profiles is ejected to radial shells 
directly outside the cusp during a merger 

6. Particles in a given radial shell outside the cusp are re- 
distributed almost uniformly with radius primarily due 
to phase mixing. This result holds for both the isolated 
spherical halo and the merging halos. The merger helps 
to redistribute an additional ~20% of orbits that have 
predominantly tangential velocity distributions in the 
original halos. 

7. Particles within the scale radius of the merging NFW 
halos experience a net decrease in total energy follow- 
ing the merger. This is a result of the overall decrease 
in the potential energy of the particles, resulting from 
the (~ 60%) increase in the mass of the cusp in the final 
remnant. 

8. The majority of particles in the merging halos expe- 
rience only a small change in their kinetic energies 
and moderate changes in their potential energies. The 
change in total energy of particles is driven predomi- 
nantly by the changes in their potential energies during 
pericenter passages. 

9. At smaller radii the majority of particles experience 
only small change in angular momentum evidenced by 
highly peaked distributions of AJ. The angular momen- 
tum of the merger is absorbed primarily by particles in 
the outermost radial shells. 

10. The largest fraction of particles lying beyond the virial 
radius of the final remnant (25%) were located beyond 
the virial radius of the original halos. The inner 2 shells 
are most robust to the ejection of particles beyond the 
virial radius. All shells between shell 4 and the virial 
radius (shell 10) contribute roughly equally to the par- 
ticles ejected beyond /?vir of the remnant. 

While this paper focuses on simulations of the merger of 
two dark matter halos with NFW potentials, most of the re- 
sults are generic to mergers of collisionless gravitating sys- 
tems and therefore applicable to dissipationless ("dry") merg- 
ers of elliptical galaxies. In particular, the absence of large 
amounts of mixing in radius over and above that expected 
from phase mixing in an isolated halo supports previous work 
that indicates that radial gradients in stellar properties (such as 
metallicity gradients}^ are can to survive intact even in major 
mergers dWhitellTQSOl: iBarnesiri996t iBovlan-Kolchin & Mai 
|2004|) . These results are also expected to be generically 
applicable to cosmological mergers. In fact soon after the 
submission of this paper we became aware of the work of 
Faltenbache r et al.l(l2006h - a study of the relaxation processes 
that operate in hierarchical mergers of dark matter halos in 
cosmological A^-body simulations. These authors indepen- 
dently arrived at the same conclusions that we do, namely that 
the primary driver for mixing and relaxation in cosmological 
mergers is tidal shocking. 

Finally, we find strong evidence that the processes that drive 
the mixing in phase space and evolution to a dynamical equi- 
librium do not occur continuously due to the time-dependent 
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potential of the system. In fact strong mixing occurs primar- 
ily following episodic injection of energy and angular mo- 
mentum into the internal energies of particles during pericen- 
ter passages. The injection of energy and angular momentum 
causes nearest particles in phase-space to separate in E, J. 

Previous studies of orbital evolution in time-dependent and 
time-independent potentials have argued that the presence of 
large fractions of chaotic orbits could be the principal driver 
of the mixing and evolution to an equilibrium distribution. 
While chaotic mixing could well be occurring, we do not find 
evidence that such chaotic mixing is driving the relaxation. 
This is probably a consequence of the fact that the timescales 
required for chaotic mixing to be relevant are much longer 
than the duration of the merger. 

Two interesting questions remain regarding "violent relax- 
ation" 1) does the mixing that occurs during collapse of an 
isolated cold gravitating system give rise to similar mixing in 
phase-space? 2) what gives rise to the universal power-law 
phase-space density profiles seen in cosmological mergers? 
We defer addressing both questions to future work. 
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